## Reproduction Script for:
## Robles, Pedro and Daniel J. Mallinson
## "Policy Learning and the Diffusion of Autonomous Vehicle Policy in the American States
## State and Local Government Review
## https://doi.org/10.1177/0160323X241262816

rm(list = ls())

library(lme4)
library(psych)

############################### Dyadic EHA Models

data <- read.csv("av_dyadic_data.csv")

## Examine which random effect structure best fits the data using guidance from 
## Field, Andy, Jeremy Miles, and Zoe Field. 2012.  Discovering Statistics Using R. SAGE. 

interceptonly <- glm(adopt_2~1, family=binomial(link="logit"), data=data)
summary(interceptonly)

policyintercept <- glmer(adopt_2~ 1|topics_2, family=binomial(link="logit"), data=data)
summary(policyintercept)

lrtest(interceptonly, policyintercept)

bothintercept <- glmer(adopt_2~ (1|topics_2) + (1|state_2), family=binomial(link="logit"), data=data)
summary(bothintercept)

lrtest(policyintercept, bothintercept) #Including random intercepts for both policy and state is the best fitting model

## Duration Dependence Choice

linear <- glm(adopt_2 ~ time, family=binomial(link = "logit"), data=data)

loglinear <- glm(adopt_2 ~ log(time), family=binomial(link = "logit"), data=data)

data$time2 <- data$time^2
data$time3 <- data$time^3

polytime <- glm(adopt_2 ~ time + time2 + time3, family=binomial(link = "logit"), data=data)

lrtest(linear, loglinear)

lrtest(linear, polytime)

#####
##Descriptive Table for Dependent and Independent Variables
D00 <- data[c("adopt_2","neighbor","ideodist", "totalpolicy_2", "legprofsquire_2", "poptotal_2","incomepcap_2","EPU_State_2","inst6017_nom_2","electronics_ig_2","transportation_ig_2","time")]

describe(D00)

## Simple Model

model1 <- logitor(adopt_2 ~ neighbor + ideodist + legprofsquire_2 + poptotal_2 + incomepcap_2 + EPU_State_2 + inst6017_nom_2 + electronics_ig_2 + transportation_ig_2 + time, data=data,robust=TRUE)
model1

cbind(Est=model1$oddsratio[,1], LL=model1$oddsratio[,1] - 1.96*model1$oddsratio[,2], UL = model1$oddsratio[,1] + 1.96*model1$oddsratio[,2])

## HLM
#Rescale data
scaled <- scale(D00[,2:12])
scaled <- as.data.frame(cbind(data$adopt_2, scaled,data$year_2))
scaled <- cbind(scaled, data$state_2, data$topics_2)
names(scaled)[c(1,13,14,15)] <- c("adopt_2", "year_2", "state_2", "topics_2")

model2 <- glmer(adopt_2 ~ neighbor + ideodist + legprofsquire_2 + poptotal_2 + incomepcap_2 + EPU_State_2 + inst6017_nom_2 + electronics_ig_2 + transportation_ig_2 + time + (1|topics_2) + (1|state_2), data=scaled, control=glmerControl(optCtrl=list("nlminbwrap", maxfun=5e10)), family=binomial(link="logit"))

ss <- getME(model2,c("theta","fixef"))
m2 <- update(model2,start=ss,control=glmerControl(optCtrl=list(maxfun=2e4)))

ss <- getME(m2,c("theta","fixef"))
m3 <- update(m2,start=ss,control=glmerControl(optCtrl=list(maxfun=2e4)))

ss <- getME(m3,c("theta","fixef"))
m4 <- update(m3,start=ss,control=glmerControl(optCtrl=list(maxfun=2e4)))

ss <- getME(m4,c("theta","fixef"))
m5 <- update(m4,start=ss,control=glmerControl(optCtrl=list(maxfun=2e4)))

summary(m5) #Final results

se <- sqrt(diag(vcov(m5)))
tab <- cbind(Est=fixef(m5), LL=fixef(m5) - 1.96*se, UL = fixef(m5) + 1.96*se)
round(exp(tab),2)

model1a <- glm(adopt_2 ~ neighbor + ideodist + totalpolicy_2 + legprofsquire_2 + poptotal_2 + incomepcap_2 + EPU_State_2 + inst6017_nom_2 + electronics_ig_2 + transportation_ig_2 + time, data=data, family=binomial(link="logit"))

lrtest(model1a, m5)


